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ABSTRACT 

We extend the previously described CMB Gibbs sampling framework to allow for exact Bayesian 
analysis of anisotropic universe models, and apply this method to the 5-year WMAP temperature 
observations. This involves adding support for non-diagonal signal covariance matrices, and imple- 
menting a general spectral parameter MCMC sampler. As a worked example we apply these techniques 
to the model recently introduced by Ackerman et al., describing for instance violations of rotational 
invariance during the inflationary epoch. After verifying the code with simulated data, we analyze the 
foreground- reduced 5-year WMAP temperature sky maps. For I < 400 and the W-band data, we find 
tentative evidence for a preferred direction pointing towards {l,b) ~ (110°, 10°) with an anisotropy 
amplitude of 5* = 0.15 ± 0.039. Similar results are obtained from the V-band data [5* = 0.10 ± 0.04; 
{I, b) = (130°, 20°)]. Further, the preferred direction is stable with respect to multipole range, seen in- 
dependently in both £ = [2, 100] and [100,400], although at lower statistical significance. We have not 
yet been able to establish a fully satisfactory explanation for the observations in terms of known sys- 
tematics, such as non-cosmological foregrounds, correlated noise or asymmetric beams, but stress that 
further study of all these issues is warranted before a cosmological interpretation can be supported. 
Subject headings: cosmic microwave background — cosmology: observations — methods: numerical 



1. INTRODUCTION 

Since the early 1990's, great advances have been made 
in the field of data analysis techniques for studying 
the cosmic microwave background (CMB). Observations 
of the CMB anisotropics, for instance those made by 
the Wilkins on Microwave Aniso t ropy Probe (WMAP ) 
experiment (jBennett et all l2003l : iHinshaw etall I2007D . 
provides the single most powerful probe in contempo- 
rary cosmology. From these, various theoretical universe 
models may be constrained, and today an effective con- 
cordance model based on the inflationary ACDM frame- 
work has been established. 

The theory of inflation was initially propo sed as a so- 
luti on to the horizon and flatness problem ()Guth et al I 
Il981i) . Additionally, it established a highly successful 
theory for the formation of primordial density pertur- 
bations, thus providing the required seeds for the large- 
scale structures (LSS), later giving rise to the tempera- 
ture anisotropics in the cosmic micr owave background 
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A firm prediction of infiation is that the observed 
universe should be nearly isotropic on large scales. 
Yet, recent theoretical studies have demonstrated that 
anisotropic inflatio n ary m od els are indeed conceivable 



lYokovama fc Sodal l2008l). Two oth er exan iples are 
thos e presented b y lAckerman "eTall (|2007r > (ACW) 
and lErickcek et al.l ()2008l ). The first model considers 
violation of rotational invariance in the early universe, 
while the second model describes the effects on the 
observed perturbation distribution due to a large-scale 
curvaton field. 

The introduction of anisotropic models poses several 
problems in terms of data analysis. The definition of a 
proper likelihood function may be non-trivial for a gen- 
eral case, although many models can be described as mul- 
tivariate Gaussians with non-diagonal covariance matri- 
ces. All models mentioned above are examples of this. 
Yet, even in these relatively simple cases, the numerical 
evaluation of the likelihood is computationally unfeasible 
due to the sheer size of the relevant covariance matrix. 

In the present paper, we extend the pr eviously de- 
scribed CMB Gibbs sampling framework (iJewell et al.l 
l2004HWandelt et al.ll2004HEriksen et al.ll2004b( l to allow 
for non-diagonal, but sparse, covariance matrices. As 
currently described in the literature, this framework al- 
lows for exact Bayesian analysis of high-resolution CMB 
data, but only under the assumption of isotropy, i.e., 
a diagonal CMB covariance matrix. This method has 
already been applied several times to the WMAP data 
(jO'Dwver et al.ll200llEriksen et al.ll2007ai rhl. l2008bD . and 
has been extended t o take into account both polarization 
Larson et al. 2007) and internal component separation 
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The question of isotropy has received considerable at- 
tention during recent years, due to unexpected signatures 
observed in the WMAP sky maps. These data appear to 
exhibit several significant and dis tinct signatures of viola- 
tion o f statistical isotropy. First. Ide Oliveira-Costa et al.l 
(|2004| ) found a striking alignment between the two largest 
harmonic modes in the temperature anisotropy sky, the 
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Fig. 1. — Covariancc elements, Ci^i and Cf_f-|_2, used in the 
construction of the ACW covariancc matrix. These are computed 
by modifying CAMB, a publicly available Boltzmann code. 

quadrupole and the octopole. Second, IVielva et al.l 
(|2004[ ) pointed out the presence of a very large cold spot 
in the southern Galactic sky, apparently incompatible 
with A CDM-based simulations. Finally, lEriksen et al.l 
(|2004aD found a significantly anisotropic distribution of 
power between two hemispheres. The tools developed in 
the present paper may be able to constrain specific mod- 
els relevant for these observations. In particular, we use 
these methods to estimate the anisotropy parameters in 
the ACW model from the 5-year WMAP temperature 
data. 

The paper is structured as follows: In Sj^l we review the 
ACW universe model, and briefly introduce the relevant 
posterior distribution. Next, we present the method in 
before we apply our tools to simulated data in 311 
In ^ we analyze the five-year WMAP temperature sky 
maps. Finally, we conclude in ^ 

2. THE ANISOTROPIC ACW UNIVERSE MODEL 

There has been a surge of interest in anisotropic uni- 
verse models since the release of the 1-year WMAP data 
in 2003, when several hints of violation of statistical 
isotropy and/or non-Gaussianity were reported. One 
such model was devised by ACW in order to study vi- 
olations of rotational invariance during the inflationary 
epoch. In this section, we briefly review this model as it 
will be used as an worked example of the general analy- 
sis framework. However, we emphasize that the methods 
described in this paper are general and suitable for any 
universe model that predicts a sparse CMB signal covari- 
ance matrix. 

ACW considered breaking of rotational invariance by 
generalizing the spectrum of primordial density pertur- 
bations P{k) to include a preferred direction, fi, as well 
as wave- number k, 

P(k) = P(fc)(l+g(fc)(k.n)2). (1) 

Here k is the unit vector along k. and g{k) is a general 
function of k. Using a combination of naturalness ar- 
guments and detailed analysis of specific models, ACW 
then argued that g{k) in most cases can be well approx- 
imated by a simple constant, g*, and presented the full 
CMB covariancc matrix corresponding to this modified 
power spectrum, 

Slra.t'm' = CeSu'Smm' + Airn,t'm'- (2) 



Here Sim,i'm' — {aimCi^im') ^^"^ CMB signal covariance 
matrix, Cg is the angular CMB power spectrum given as 

Cii = j dkk^P{k)e^{k) (3) 

where &e{k) is the transfer function. The term A{„i.i'm' 
is then defined as 

Aem,i'm' = gMTn,i'm' dkP P{k)Qi{k)Qi, (k) (4) 

In this expression, £,em.,i'm' are geometric coefficients (see 
ACW for explicit details). The ^ coefficients couple £ to 
f = {£,£± 2} and m to m' = {m, m ± 1, m ± 2}. AU 
other elements are zero. 

The coupling to standard cosmological parameters en- 
ter only through Ci^f, which is a straightforward gen- 
eralization of the angular CMB power spectrum. In 
this paper, we assume that the cosmological parame- 
ters are known, and only the anisotropy parameters, 
g^ and n. are unknown. We therefore compute Ci^i' 
once, using a very slightly modified version of CAMB 
(jLewis et al.l |2000[ ) that outputs Cg,g+2 in addition to 
Cf^f, and adopt this matrix as a prior. We adopt the best- 
fit ' ACm£jnodel_detemm^ from the 5-year WMAP 
data (jKomatsu et al.|[2008[ ). and the corresponding Ci^i 
and Ct,i+2 elements are plotted in Figure [1] Joint es- 
timation of cosmological parameters and the anisotropy 
parameters will be considered in a future publication. 

In Figure [2] we show one realization drawn from a 
Gaussian distribution with zero mean and Sim/'m as 
covariancc matrix, with g^ = 0.9999 (middle row) and 
= —0.9999 (bottom row) and a preferred direction of 
(/,6) = (0°,90°). The isotropic signal is depicted in the 
top row. 

The anisotropic contribution alone consists of correla- 
tions with the underlying isotropic signal stretched along 
the plane normal to the preferred direction. The sign of 
determines whether the anisotropic contribution is to 
be added or subtracted from the isotropic signal. If the 
anisotropic signal is added, then the spots are stretched 
along the plane normal to the preferred direction. How- 
ever, if the anisotropic signal is subtracted (g* < 0), then 
the spots are effectively squeezed along the plane normal 
to the preferred direction, corresponding to stretching 
parallel to the preferred direction. 

2.1. The As- g^: degeneracy 

From equation ([2]) and the definition of ^ (see ACW) 
it is clear that A contributes also to the diagonal of the 
signal covariance matrix, and therefore affects the total 
angular power spectrum, not only the correlations among 
aim's. This introduces a strong degeneracy between 
and the amplitude of the power spectrum of scalar per- 
turbations, As or (Tg. Unless one attempts to estimate 
the standard ACDM parameters jointly with the new 
anisotropy parameters, one must therefore ensure that a 
given choice of does not significantly affect the overall 
power spectrum, but only the anisotropic contribution. 

The diagonal part of A. for which the integral over the 
transfer functions equals Ci, is 

Afm,fm = g*Ce£,em,em- (5) 
The net extra power due to A is therefore 

9 Ci ^ 

-Df ^^Tfry X! ^i^n/m- (6) 
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Fig. 2. — Temperature maps showing isotropic fluctuations (top row), while the two lower rows depict anisotropic contributions with 
Qt = 0.9999 (middle row) and g» = —0.9999 (bottom row). The maps in the left column are presented in Mollweide projection, while the 
right row is Cartesian. The anisotropy direction was chosen to be {l,b) = (0°,90°). Note the subtle tendency for stripes along the equator 
for the positive g», and perpendicular to the equator for negative g*. 



This may be greatly simplified by considering the de- 
tailed form of ^lm,lra-, 



Averaging this expression over m, one finds that 
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such that Di = ig^Cf. We therefore redefine the total 
signal covariance matrix to read 



1 + 5*/3 



(10) 
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Fig. 3. — Power spectra of simulated anisotropic sky maps with 
Qt = 3, with (green) and without (black) rescaling. Red curve 
shows the power spectrum for an isotropic simulation with g, = 0. 

With this definition, is a direct measure of the 
anisotropic component of S, and docs not directly de- 
pend on the power spectrum Ci. 

The effect of on the power spectrum is demonstrated 
in Figure [3l where we plot the power spectra of a simu- 
lated anisotropic map with 17* = 3, with and without the 
above rescaling. Unless proper rescaling is performed, or 
some equivalent parametrization introduced, it is clear 
that the strongest constraints on will come from the 
observed power spectrum, rather than the correlations 
among aim's. 

2.2. Posterior analysis and priors 

The goal is now to estimate 5, and n from observed 
CMB maps, by computing the posterior distribution 
P((7*,n|d), d denoting the data. Because we assume 
that both the noise and CMB sky signal are Gaussian 
(but anisotropic) random fields, this distribution reads, 
by Bayes' theorem, 

P(g„n|d)«£(5„n)P(5„n), (11) 

where £(5,, n) = F(d|5,, n) is the likehhood 

/:(.g„fi)cx y=- (12) 

V 1^1 

and P((7», fi) is a prior. Equatio n (IT^ can be eva luated in 
C(^pix) operations, as shown bv lOh et al1()1999D . In this 
expression, C is the signal-plus-noise covariance matrix. 
In principle, we could now simply map this distribution 
over a three-dimensional grid, and our task would be 
completed. However, except for the special case of a 
data set with uniform noise and full-sky coverage, this 
is in practice impossible because C is a dense matrix, 
and inversion and matrix determinant therefore scales as 
^(^pix)' ^pix being the number of pixels. For current 
and future data sets, one expects A^pix ~ 10^ or more. 

Fortunately, there is one specific feature of the ACW 
model that does make an exact analysis possible: Al- 
though the full-sky CMB covariance matrix is non- 
diagonal, it is not dense. Rather, it has a well-defined 
shape in harmonic space {£ is coupled to £' = {1,1 ± 2} 
and m to m' = {to, to ± 1, to ± 2}) that allows for cheap 
matrix storage and fast Cholesky decomposition. This, 
combined with the development of the standard diago- 
nal CMB Gibbs sampler mentioned in the introduction 
(iJewell et all l2004t IWandelt et a L '2004'; lEriksen et aP 
l2004b[ ). allows us to perform a full proper analysis, as 
explained in the next section. 
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Fig. 4. — Marginal likelihood functions for jC(g,) (top) and C{ri) 
(bottom) for a simulated data set with uniform noise and full-sky 
coverage, shown in logarithmic units. The input values of 3» = 
0.8 and (/,fe) = (57°, 33°) are accurately reproduced. Notice the 
shallow local maximum at g« ~ —0.5 and the secondary peaks in 
the marginal direction map. 

Before describing this method, we consider first the 
special case of data having uniform noise and full-sky 
coverage, which is useful to illustrate the approach, and 
highlight some particular issues. For this particular case, 
the full data covariance matrix, expressed in spherical 
harmonic space, has the same sparse filling pattern as the 
ACW covariance matrix, and direct evaluation is there- 
fore p ossible using sparse matrix techniques (^e.g.. [Pavisl 



We simulated a single CMB realization from the ACW 
model, adopting a high anisotropy amplitude of — 0.8 
and a preferred direction (in Galactic longitude and lat- 
itude) of {l,b) = (57°, 33°), then convolved this real- 
ization with a 90' FWHM Gaussian beam, and pro- 
jected it onto a HEALPix^ grid with resolution param- 
eter A^sidc = 128. Finally, uniform, Gaussian noise with 
lO^K RMS was added to each pixel. This simulation 
was then analyzed by computing the raw likelihood over 
a three-dimensional grid, and finally marginalized like- 
lihoods were produced by numerical integration. The 
results from this exercise is shown in Figure |4l 

As expected, the likelihood peaks close to the input 
values. However, there is also a second local maximum 

at 5* 0.5 with a direction of (I, b) ~ (45°, -50°), 90° 

with respect to the main axis. This maximum becomes 
visible only for large negative values of . The existence 
of this maximum becomes intuitive when considering fig- 
ure [21 Flipping the sign of and rotating the preferred 
axis by 90° leads to stripes in the same direction as the 
original parameters. 

This is not a significant issue for a direct evaluation 
method, since the local maximum has a very small am- 
plitude. (Note that the marginal likelihoods in Figure 

http://healpix.jpl.nasa.gov 
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|4]are shown in logarithmic units.) However, for MCMC 
methods it can cause problem in terms of burn-in: As 
explained in the next section, our method is based on 
the well-known MCMC and Gibbs sampling algorithms, 
and these essentially correspond to performing a random 
walk on the likelihood surface. Further, each chain is 
initialized randomly on the sphere. It is therefore a sig- 
nificant chance that a number of chains may get trapped 
in a local maximum, and thereby bias the final posterior. 
To avoid this, we impose a uniform prior of > —0.2 
in this paper, and a uniform prior on the sphere for n. 
If the final posteriors from the actually WMAP analysis 
happened to peak close to = —0.2 we would have to 
re-consider this choice more carefully, but as we shall see, 
this is not the case. 

3. METHOD 

We now discuss the method for mapping out the de- 
sired posterior. This method is a very slight general- 
ization of the pre viously descr i bed C M B Gibbs sam- 
pler d evelo ped bv iJewe U et all (120041) ■ IWandelt et all 
()2004f) and Eriksen et all (12004^ ). which was originally 
intended for power spectrum estimation. The underly- 
ing Gibbs sampler implementation used for this work 
is the code ca l led "Comman der", described in detail by 
lEriksen et all (I2004E [2008al) . 

3.1. Review of the CMB Gibbs sampler 

We first review the CMB Gibbs sampler as previously 
described in literature. In any Bayesian analysis, a main 
goal is the posterior distribution P{9\d), where is a set 
of parameters connected to some model and d are the 
observed data. For high-dimensional spaces, brute-force 
evaluations of the posterior are computationally unfeasi- 
ble, and one usually resorts to Monte Carlo Markov chain 
(MCMC) methods. 

3.1.1. Notation and data model 

We begin by defining a parametric model for the CMB 
observations. Given our current understanding of the 
CMB sky, the observed data may be accurately modelled 
as a sum of a CMB anisotropy term and a noise term. 



d = As + n. 



(13) 



Here d represents the observed data, A denotes 
convolution by an instrumental beam, s(0, 0) = 
m 'z^) ^^'^ CMB sky signal represented 
in either harmonic or real space, and n is instrumental 
noise. 

Further, it is a good approximation to assume both the 
CMB and noise to be zero mean Gaussian distributed 
variates, with covariance matrices S and N, respectively. 
In harmonic space, the signal covariance matrix is de- 
fined by Sem,i'm' = {aimae'rn')^ which may or may not 
be diagonal. The connection to cosmological parameters 
is made through this covariance matrix. Finally, for 
experiments such as WMAP, the noise is often assumed 
uncorrelated between pixels, N^j = for pixels i and 

j, and noise RMS equals to ct^. 

Our goal is now to compute the full joint poste- 
rior P{9\d), which, as already mentioned, is given by 
Pi9\d) cx F(d|e')P(6') = C{e)P{e), where C{e) is the 
likelihood, and P{9) is a prior. For a Gaussian data 



model, the likelihood is 
C{9) oc 



(14) 



3.1.2. Posterior mapping by Gibbs sampling 



When working with real-world CMB data, there are 
a number of issues that complicate the analysis. Two 
important examples are anisotropic noise and Galactic 
foregrounds. First, because of the scanning motion of a 
CMB satellite, the pixels in a given data set are observed 
by unequal amounts of time. This implies that the effec- 
tive noise is a function of position on the sky. Second, 
large regions of the sky are obscured by Galactic fore- 
grounds (e.g., synchrotron, free- free and dust emission), 
and these regions must be rejected from the analysis by 
masking. 

Because of these issues, the total data covariance ma- 
trix S + N is dense in both pixel and harmonic space. 
As a result, it is computationally difficult to evaluate 
the likelihood in Equation (fT4|) . since the computational 
cost of matrix inversion and determinant evaluation scale 
as C'(A''p;^). Fortunately, this problem has already been 
solved for the CMB context, through the development of 
the CMB Gibbs sampler. 

The idea behind the CMB Gibbs sampler is to estimate 
the CMB sky, s, together with the covariance parameters, 
by computing P{9^ s|d), and then subsequently marginal- 
ize over s. Specifically, the algorithm is the following: 
First choose any initial guess, {9,sf. Then alternately 
sample from each of the conditional distributions, 

9'+^*- P{9\s\d) (15) 

s*+i ^P(s|e''+i,d). (16) 

The theory of Gibbs sampling then guarantees that the 
joint samples (0, s)* will, after some burn-in period, be 
drawn from the desired joint distribution. The remaining 
step is then simply to formulate sampling algorithms for 
each of the two conditionals, P{9\s, d) and P{s\9, d). 

We first consider P(s|6',d). This may. under the as- 
sumption of Gaussianity, be written as 

P(,s\9,d) = P{d\s,9)P{s\9) (17) 



,-i{d-s)^N-i(d-s)g-isS- 
^-i(s-^S)^(S-i+N-i)(s-S) 



(18) 

(19) 

where we have defined the Wiener filtered map, s = 
(S-i + N-i)-iN-id. Thus, P(s|6l, d) is a Gaussian dis- 
tribution with mean s and covariance (S~^ -I- N~^)~^. 

Sampling from this distribution is straightforward, but 
implementationally somewhat involved: Draw two Gaus- 
sian random maps, 770 and 771, with zero mean and unit 
variance, and solve the following equation for s, 

(S-i +N-i)s = N-M-hL-^r/o + N-3, (20) 

where L is the Cholcsky decomposition of S = LL-^. 
By multiplying both sides of this equation with (S^^ + 
N^^)^^, one immediately sees that (s) = s, and a 
few more computations show that ((s — s)(s — s)^) = 
(S^^ + N~^)~^. as required. 

For improved numerical stability, this linear system is 
in practice rewritten into the following form, 

-L'^N"^ (21) 



(1 + L^N"iL)(L"is) = L'^N^M + ?/o - 
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which is first solved for x — L^^s by conjugate gra- 
dients, and then for s = Lx. For fur t her im plemen- 
tational details, see, e.g., lEriksen et al.l (|2008al ). Note, 
however, that in previous papers equation (|2ip was al- 
ways written with symmetric signal covariance square 
roots, S2 = (82)"^. The current form is based on the 
Cholesky decomposition, which is computationally con- 
siderably cheaper than the symmetric form, especially 
for sparse matrices. 

Finally, we need a sampling algorithm for P(0|s,d). 
In previous publications, the main emphasis has been on 
covariance matrices parametrized by the angular CMB 
power spectrum, Cim/'m' = CiSu'Smm' ■ In this case, 
P(Cf|s, d) reduces to a simple inverse Gamma distri- 
bution, for which there is a simple textbook sampling 
algorithm available. As the details of this specific algo- 
rithm is of little use for the application presented here, 
we refer the interested reader to earlier pape r s for full 
details on this pr o cedur e, e.g.. IWandelt et all ()2004l ) or 
lEriksen fc Wehi^ (I2OO8I ). 

3.2. Gibbs sampling with non-diagonal covariances 

We now describe the two modifications to the CMB 
Gibbs sampler that allow us to analyze models with non- 
diagonal covariances. This involves adding support for 
non-diagonal covariance matrices for P{s\6, d) and imple- 
menting a more general sampling algorithm for P(^|s, d). 

3.2.1. Sampling from P{s\0, d) 

We first consider sampling of sky maps, s, given a set 
of cosmological parameters, 6, and the associated co- 
variance matrix S{9). Formally, the sampling algorithm 
for P{s\9,d) is identical to that given by equation ([2T|) . 
However, in this case S is a non-diagonal matrix and the 
computational complexity is therefore greatly increased. 
Only special cases can be considered, for instance models 
that predict a sparse covariance matrix. This is the case 
for the ACW model. 

For general dense anisotropic covariance matrices, the 
memory requirements scale as 0{£'^^^^), effectively ren- 
dering studies of anisotropic models where £max ^ 100 
impossible. However, working only with sparse matrices, 
the memory consumption scales as ©(-^max): enabhng cal- 
culations of covariance matrices with ^max well into the 
Planck regime (^max ~ 2500). 

To be able to handle sparse matrices effi ciently, we 
have ported the LDL library of iDavisI (|2005l ) to Fortran 
90, and incorporated this into Commander. This library 
stores sparse matrices in a packed format, and supports 
fast Cholesky decomposition. Our F90 version of LDL 
may be obtained by sending an email to the authors, and 
will be released publicly at a later time. 

In the present paper, we are primarily concerned with 
the ACW model, and the corresponding covariance ma- 
trix exhibits correlations between £ and £' = {£, £±2} and 
between m and m' = {m, to± 1, m± 2}. Thus, the num- 
ber of elements up to fmax is 0{15£^^^). For example, 
for ^niax = 300 the memory requirements are ~ 14Mb 
with double precision complex numbers. Since the co- 
variance matrix is very sparse, the CPU time required 
for Cholesky decomposition is nearly linear in £,^„ax- 

We define three different harmonic space limits in our 
code, namely €max, £iow and ^high- The former denotes 
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Fig. 5. — Evolution of Gibbs chains mapping the posterior of 
a simulated data sot. Note how chains trapped in the local maxi- 
mum at negative anisotropy amplitude eventually converges to the 
positive maximum. 

the maximum multipole moment of the full spherical har- 
monics composition used in the analysis, while the latter 
two denotes the range in which the anisotropic covariance 
matrix is used. In addition, we remove the monopole and 
dipole from the analysis. Thus, the total covariance ma- 
trix reads 





mm' 



i+g,/3 



£,£'<! 
2<£,£' <£io^ 

-^low — ^ — -^high 
-^high *\ ^ -^max 

(22) 

The reason for defining ^min and ^max as free parame- 
ters is that it may be useful to study the dependence of 
6* on a particular ^-range. On the other hand, this im- 
plies that the model implemented in this paper is only 
an approximation to the full ACW model, for which the 
correlations extend over all fs, and is only exact when 

^high — ^max- 

With the sparse matrix operations implemented, the 
algorithm is precisely the same as for the diagonal case, 
and both rely on the s olution of a linear sys tem by Con- 
jugate Gradients fCG: lEriksen et ani2004b[ ). In order to 
achieve an acceptable CG convergence rate, it is there- 
fore necessary to establish a good preconditioner. How- 
ever, as long as the off-diagonal elements remain small, 
the standard diagonal covariance matrix preconditioner 
performs reasonably even for the off-diagonal case. For 
the present paper, we therefore adop t the sa me precondi- 
tioner as described bv lEriksen et al.l ()2004bf ). which con- 
sists of the directly inverted full matrix evaluated up to 
some ^precond, and then a strictly diagonal matrix from 

^prccond ~t- 1 tO ^max- 

The number of CG iterations per map making step is 
typically 70 for a WMAP-type run, and with a total CPU 
time per iteration of about 15 seconds, the total cost 
for a single sample is ~ 20 CPU minutes. The average 
CPU time required to set up and perform a Cholesky 
decomposition of the corresponding covariance matrix for 
^max = 512, flow = 2 and £high = 300 is ~ 20 seconds. 

3.2.2. Sampling from P (6 \s,d) 

Finally, we have to formulate a sampling algorithm for 
P{9\s,d). Recall that for the diagonal power spectrum 
case, this step is typically performed b y a standard in- 
verse Gam ma distrib ution sampler fe.g.. lGupta fc Nagail 
HOOO: Erik sen fc Weh us 2008i). For the general case con- 
sidered here, we adopt a standard Metropolis MCMC 
sampler fe.g.. lLiul[200l[) . 

First, note that P{6\s, d) = P{6\s); if we already know 
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Fig. 6. — Posterior distributions for simulated maps with a significant anisotropic amplitude g* = 
gt = 0.0 (right). Note how the anisotropic input parameters were successfully reproduced. 
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the CMB sky perfectly, no additional data can possibly 
tell us anything more about the anisotropy parameters 
9. Second, although the CMB sky is now manifestly 
anisotropic, we still assume that it is Gaussian, and the 
target distribution therefore reads 



p{e\s) cx 



(23) 



For sparse matrices, this may be directly evaluated by 
first computing the Cholesky decomposition of S = LL , 
and then, on the one hand, solve for x = Ls, and on the 
other hand, compute |S| = |Lp. 

We adopt a simple symmetric proposal rule for the 
Metropolis sampler, and the acceptance probability 
therefore simply reads 



P 



P{9P 



(24) 



where 9^ is the proposed sample and 0* is the current 
sample of the MCMC chain. Specifically, we adopt a 
Gaussian proposal density for and a uniform proposal 
over a disk for n, centered on the current state. The pro- 
posal density is typically tuned by producing a short test 



chain before the main run, such that the final observed 
acceptance rate lies between 0.2 and 0.7. 

Finally, because the computational cost is much lower 
for this step than for P{s\9,d), we produce several 9 
samples per main Gibbs iteration, to improve the con- 
vergence properties of the chain. This essentially cor- 
responds to per forming a partial Rao-BlackwcUization 
(|Chu et al.ll2005D . A typical number of MCMC samples 
per main Gibbs iteration is 30. 

4. APPLICATIONS TO SIMULATED DATA 

We now apply the methods described above to sim- 
ulated data, both in order to validate the code and to 
build up intuition about the target distribution. Note 
that the discussion from now on specializes exclusively to 
the ACW model, and it is possible that other technical 
issues than those described here may arise when consid- 
ering other models. Burn-in, mixing and convergence are 
issues that must be considered on a case-to-case basis. 

4.1. Simulations 

To test our implementation and study the behavior 
of the algorithm in general, we simulate a few different 
maps from the ACW model, and analyze these maps with 
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Fig. 7. — Posterior distributions for a simulated WMAP data 
set, using the V-band beam, V band RMS noise and the KQ85 sky 
cut. 

our modified Gibbs sampler. The CMB component of 
these maps is made by generating a random vector, 77, 
of Gaussian uniform variates with zero mean and unit 
variance, and then computing s = L77. This realization is 
then convolved with a beam function and the HEALPix 
pixel window, before it is projected on a HEALPix grid. 
Finally, Gaussian noise is added to each pixel. 

The first two simulations have a resolution of A^sidc = 
128, €,„ax = 256, iiovj = 2, ihigh = 200 and a Gaussian 
beam of 90' FWHM. The noise RMS is 10 fiK uniformly 
over the full sky. The first of the two simulations has an 
anisotropy amplitude of 5, = 0.8 and a preferred direc- 
tion towards {l,b) = (57°, 33°), and the other = 0. 
These two simulations are primarily used to compare 
the Gibbs sampler with brute-force likelihood evaluation, 
which is only possible for uniform noise and full-sky cov- 
erage. 

Second, we generate a full WMAP5 like simulation 
based on the VI differencing assembly (DA), with = 
0.8, iV,ide = 512, Cax = 600, ^low = 2, fhigh = 300 and 
beam and noise properties appropriate for t he VI DA*^. 
In th is case, we also apply the KQ85 sky cut (jGold et al.l 
[2008I) . which removes 18% of the sky. This simulation 
is used to verify that correct results arc obtained for re- 
alistic WMAP data, including anisotropic (but uncorre- 
lated) noise and a sky cut. 

4.2. Burn-in and convergence 

We first consider the issue of burn-in and convergence, 
and analyze the simulation with 5, = 0.8, uniform noise 
and full-sky coverage. In Figure [5] we show the first 6000 
(7* samples produced by each of 14 chains. First, notice 
that the chains immediately divide into two classes, one 
which converges quickly towards 5, ~ 0.8 and one which 

http://lambda.gsfc.nasa.gov 



hovers near the lower prior of 5, — —0.2. This is due 
to the fact that the chains are initialized randomly on 
the sphere, and those that happen to start close to the 
non-physical local maximum (sec Section 12. 2[) get tem- 
porarily trapped in this local maximum. However, as 
the chains explores the likelihood surface, they are able 
to converge into the right regime, and find the correct 
value. In this case, all chains have reached the equilib- 
rium state after 1800 iterations. The pre-burn-in samples 
must be rejected from the further analysis. For now, we 
inspect each chain individually, to make sure that they 
have all reached the common state. 

Note that there is a fundamental difference between 
low and high signal-to-noise cases in this respect: If is 
low, the chains may jump between local maxima, while 
if (7* is high, some chains typically start out in the global 
maximum and stay there, while others start in the lo- 
cal maximum, and eventually converge into the right 
regime. Which situation is relevant for a particular data 
set must be considered on a case-by-case basis, by check- 
ing whether the chains jump between states, or if they 
stay in one place. It is also advisable to run many chains 
in parallel, randomly initialized over the full sphere, to 
understand how many local maxima the distribution has. 

Second, once the chains have burned in, we must 
also ensure that they collectively have converged to 
the full posterior. One possible m easure for this is 
the Gelman-Rubin R st atistic (e.g., iGelman fc RubinI 
ll992HEriksen et al.ll2006h . which compares the variances 
within a single chain with the variance between chains. 
If the chains have converged properly, R should be close 
to unity. Typically, one recommends that R should be 
less than 1.1 or 1.2. For the chains shown in Figure [5l we 
find that R = 1.01 after rejecting the first 2000 burn-in 
samples, indicating very good convergence. Considering 
further subsets of these samples, we find that 5000 sam- 
ples is sufficient to achieve i? < 1.1 and 20 000 samples 
for R < 1.02. In the analyses presented later, we always 
have more than 20 000 samples. 

4.3. Validation 

We now analyze the two simulations with uniform noise 
and full-sky coverage, having 5, = and = 0.8, respec- 
tively. In addition to running the Gibbs sampler on these 
simulations, we also compute the full three-dimensional 
likelihood function over a grid in (g» , n) , and numerically 
integrate to produce brute-force marginal posteriors. 

The resulting distributions are shown in Figure [6l the 
left column shows the g^, = 0.8 case, and the right column 
shows the = case. We see, as expected, that the 
two methods produce identical results, up to sampling 
uncertainty and grid resolution. Note that this holds 
both for high and low anisotropy amplitudes, indicating 
that the method is robust in all regimes. 

Next, we see that when the amplitude is large, there is 
only one visible preferred direction in the direction pos- 
terior; the secondary direction is too shallow to be seen. 
On the other hand, there are two "preferred" directions 
in the g^ = case. However, the span in likelihood over 
the full sphere is in this case only a factor of two between 
the least and most preferred directions, which essentially 
indicates a uniform distribution. 

In Figure [7] we show similar plots for the WMAP simu- 
lation with uncorrelatcd noise, based on the VI differenc- 
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Fig. 8. — Estimated uncertainty in 3, as a function of ^high 
(black dots) and a best-fit power law function (red line) for cosmic 

m§^mS(Mi^f^cM''g* = 0.8. in this case it is not possible 
to evaluate the likelihood directly, since the noise is in- 
homogeneous and there is a sky cut. Still, we see that 
correct results are obtained. This concludes the verifi- 
cation of both the method and our implementation, and 
we are now ready to analyze the five-year WMAP tem- 
perature sky maps. 

4.4. Forecasts for cosmic variance limited data 

Before turning to the analysis of the actual WMAP 
data, we compute the uncertainty in 5* as a function of 
^high for full-sky noiseless data. (The lower limit is al- 
wa ys kept at £iow = 2.) Th i s topi c was also considered 
bv iPuUen fc KamionkowskH (|2007f ). who presented both 
a more general formalism and forecasts for specific ex- 
periments. Note, however, that our parametrization is 
slightly different from theirs, as we introduce a rescaling 
of the covariance matrix to eliminate the power spectrum 
degeneracy fSection l2.1|) . 

We carry out this analysis by simulating anisotropic 
ACW maps with = and different ^high, and analyze 
these with the brute-force evaluation approach described 
above. No noise or beam effects are included. For each 
case, we marginalize over fi to obtain P(g, |d), and com- 
pute the standard deviation, cr(g*), from this distribu- 
tion. 

Figure [5] shows cr((7*) as a function of ^high- From this 
figure, we sec that cr(g*) is very close to a power law in 
£viiffh, in good agreemen t with the arguments given by 
iPuUen fc Kamionkowskl (|2007t ). The best-fit power law 
function is 

^(4igh;5*)- 0.025 (^^) (25) 

and this can be used to produce rough forecasts for vari- 
ous experiments. For instance, in this paper we conserva- 
tively adopt ^high = 400 for the WMAP analysis, to avoid 
possibly complicating high-.^ issues such as point source 
confusion and noise mis-estimation. In that case, we ex- 
pect an uncertainty of <j{g*) = 0.025, before taking into 
account noise and sky cut. This is in excellent agreement 
with t he cr((?*) = 0.024 result of lPuUen fc KamionkowskH 
(|2007f ). derived with slightly different data and model 
assumptions and a completely different approach. 

5. APPLICATION TO THE FIVE-YEAR WMAP DATA 

We now analyze the five-year WMAP data, and present 
the full marginal P{gt\d) and P(n|d) posteriors for var- 
ious data cuts. 



5.1. Data 

In this paper, we consider the five-year WMAP tem- 
perature sky maps (jHinshaw et al.l l2008f ) , and analyze 
the V- and W-bands (61 and 94 GHZ), which are believed 
to be the cleanest WMAP bands in terms of residual fore- 
grounds. We adopt the template-corrected, foreground 
reduced maps recommended by the WMAP team for cos- 
mologi cal analysis, and i mpose both the KQ75 and KQ85 
masks (|Gold et al.ll2008i ). which remove 28% and 18% of 
the sky, respectively. Point source cuts are imposed in 
both cases. 

We mainly analyze the data frcquency-by-frequency, 
and consider the combinations V1-I-V2 and Wl through 
W4. In addition, we compute the posteriors for VI and 
V2 separately. The noise RMS patterns and beam pro- 
files are taken into account for each DA individually. The 
noise is assumed uncorrclatcd between pixels and bands. 
For detai ls on joint Gibbs ana lysis of multi-frequency 
data, see lEriksen et al.l (|2004br ). All data used in this 
analysis are available from LAMBDA. 

The angular resolutions of the V- and W-bands are 
0.35° and 0.22°, respectively, and the sky maps are pix- 
elized at a HEALPix resolution of iVgido = 512 with 7' 
pixels. We therefore adopt a harmonic space cutoff of 
imax = 700 and 800 for the two data sets, probing deeply 
into the noise dominated regime. However, we never con- 
sider multipoles at £ > 400 for the anisotropic part of the 
signal covariance matrix, in order to minimize the chance 
of systematic effects such as residual point source contri- 
butions, beam uncertainties or noise mis-estimation to 
affect our results. See Table [T] for a list of the specific 
^-ranges considered. 

We also note that the maps stud ied here are cleaned us- 
ing external templates (jGold et al. 2008 ) , which must be 
considered a fairly rough approach to foreground clean- 
ing. A better approach is to use the joint for eground 
and CMB Gibbs sampler (jEriksen et al.l l2008a[) . which 
provides the user with a CMB map marginalized over 
very general foreground models. This work is currently 
underway for the five-year WMAP data, and the results 
will be reported elsewhere (Dickinson et al., in prepara- 
tion). However, as an explicit foreground test we also 
analyze the raw V-band data, from which no foreground 
templates have been subtracted, and find very consistent 
results. 

5.2. Results 

We now present the marginal posteriors for the ACW 
model obtained from the five-year WMAP temperature 
sky maps, as computed with the method described in fJ31 
First, in the top row of Figure |H| we show the marginal 
anisotropy amplitude posterior, P{g^,\d), for V- (left col- 
umn) and W-band (right column), and in the three bot- 
tom rows we show the preferred direction posteriors, 
P(n|d). In Table [1] the full set of results are summa- 
rized quantitatively. 

First, we sec that there is an apparently clear detection 
of 7^ when considering the full range of multipoles, 
e = [2, 400]. The W-band posterior has = 0.15 ± 0.04, 
nominally corresponding to a 3.8(t detection, and the V- 
band posterior has = 0.10 ±0.04, internally consistent 
with W-band at ^ Icr. Second, the direction posteriors 
indicate a clearly preferred direction pointing towards 
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Fig. 9. — Marginal ACW posteriors obtained from the V- (left) and W-band (right) WMAP temperature sky maps. Top row shows 
P(c;, |d) and bottom three rows show P(n|d) for three different ^-ranges. Note the common preferred axis in both £ = [2, 100] and [100, 400] . 



(/,6) = (110°,10°). 

Further, this same direction is observed in both I = 
[2 - 100] and £ = [100,400], indicating that the struc- 
ture is present over a large range of angular scales. The 
results are also stable with respect to sky cut, as the 
same pattern is seen with the KQ75 sky mask as with 
the KQ85 cut, removing an additional 10% of the sky. 

5.3. Sensitivity to systematics 

Given the nominally strong results found in the pre- 
vious section, it is imperative to search for possible sys- 
tematic effects that might explain the observations. In 
particularly, three major sources of uncertainty should 
be considered in detail, namely non-cosmological fore- 
grounds, correlated noise and asymmetric beams. 

First, residual Galactic foregrounds do not a priori ap- 



TABLE 1 

Summary of marginal posteriors from WMAP5 



Band £ range Mask 



Amplitude g« 



Direction [l, b) 



V 
V 
V 
V 

V-raw 

VI 

V2 

W 

W 

W 



2 - 400 KQ85 
100 - 400 KQ85 



2 - 100 
2 - 400 
2 - 400 
2 - 400 
2 - 400 
2 - 400 



KQ85 
KQ75 
KQ85 
KQ85 
KQ85 
KQ85 



100 - 400 KQ85 
2 - 100 KQ85 



0.10 ±0.04 
0.09[0.084, 0.148] 
-0.07[-0.156, 0.480] 
0.10[-0.100, 0.158] 
0.11 ± 0.036 
0.12 ±0.041 
0.08 ± 0.044 
0.15 ± 0.039 
0.14[-0.097, 0.236] 
0.14[-0.162, 0.470] 



(130°, 10° 
(130°, 10° 
(130°, 15° 
(130°, 10° 
(130°, 10° 
(130°, 10° 
(130°, 10° 
(110°, 10° 
(110°, 10° 
(125°, 20° 



Note. — In cases with no significant detection, the values for 
indicate the maximum posterior value and 95% confidence regions. 
Otherwise, they indicate posterior mean and standard deviation. 
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pear as a particularly promising candidate, given that 
the results are robust with respect to both frequency and 
sky cut, and the preferred axis does not point towards 
any natural Galactic direction. Second, in figure [10] we 
show the posteriors obtained from the raw V-band 5 year 
sky maps. It should be clear that galactic foregrounds 
have little impact on these results. Third, extragalactic 
point sources do also not seem as a likely candidate, be- 
cause the signature is seen both at low and high £'s, and 
we never consider multipoles above I > 400, precisely to 
avoid this type of concerns. 

The effect of correlated noise is harder to rule out. 
On the one hand, returning to the simulated anisotropic 
CMB realization in Figure [21 we see that the main sig- 
nature of the ACW model is smoothed structures along 
the plane normal to the preferred direction, and essen- 
tially no modifications along the preferred direction. On 
the other hand, the main signature of correlated noise is 
striping along the scan direction. 

Next, the ecliptic north pole has Galactic coordinates 
{l,b) = (96°, 30°), which is ~ 24° (32°) away from the 
preferred direction for W-band (V-band) found in ^15. 21 
The probability of obtaining such a close alignment by 
chance is ~ 10%, (^ 16%) which is low enough for cor- 
related noise to be considered relevant for this particular 
case. 

To study the magnitude of this effect on g*, we analyze 
realistic V and W-band simulations with correlated noise. 
These noise simulations were produced and published by 
the WMAP team in their 1-year data release. To mimic 
realistic 5-year simulations, we coadd five independent 
realizations for each differencing assembly. These noise 
realizations are then added to an isotropic CMB sky re- 
alization, and the sum is analyzed using the same proce- 
dure as in ^5.2^ We also analyze two single 1-year W4- 
band simulations, which serve as a worst-case scenario, 
as the knee frequency of this band is by far the highest 
of any WMAP DA, and the overall noise level is higher 
by a factor of ~ 4.5 than the full 5-year W-band data. 

The results from these analyses arc shown in Figure 
[TTl The left and middle columns show the simulated 5- 
year posteriors for V and W-band, respectively, and the 
right column shows the 1-year W4 posterior. The top 
row shows P(g*|d), and the bottom panel shows P(n|d). 

First, note that with realistic 5-year noise no detection 
is made in either the V- or W-band. Further, the peak 
sky position is different in the V- and W-bands, and both 
have a very low significance. It therefore seem unlikely 
that correlated noise can explain the results found in fJOj 

Still, caution is warranted, as the 1-year W4 posteri- 
ors do exhibit traces of anisotropic contributions, with 
a peak amplitude larger than the observed in the ac- 
tual 5-year data. Yet, the match with the structures 
observed in the real data is less than striking. First, the 
correlated noise simulations show two independent peaks 
in the directional posterior, while the WMAP data show 
one. Second, a detailed study of the joint posteriors for 
the simulations show that the peak along the ecliptic 
north pole corresponds to the negative peak in P(f/*|d), 
while the WMAP data has a positive 5* along its pre- 
ferred direction. Third, the preferred axis found in the 
WMAP data are further away from the ecliptic pole than 
the corresponding peak in the simulation posteriors. 

Finally, in order to make a complete analysis, we 



should also consider the impact of asymmetric beams. 
Ideally, one would prefer to address this issue in the same 
manner as correlated noise, by analyzing simulated CMB 
realizations with asymmetric beams. Unfortunately, we 
do not have access to such simulations at this time, and 
it is difficult to do a rigorous analysis. However, there are 
some arguments against the asymmetric beams hypoth- 
esis. First, the effect is observed both at low and high 
fs, with very consistent positions. Second, the observed 
preferred axis is ~25-30 degrees away from the ecliptic 
pole, and the posterior ratio of the ecliptic poles to the 
maximum posterior is low. Finally, similar signatures are 
observed in both the V and W bands, and in VI and V2, 
which all have slightly different beam patterns. 

Nevertheless, at this point it would unwise to make 
strong claims concerning a possible cosmological inter- 
pretation of the signature found in ^15.21 Proper analysis 
of fully realistic 5-year WMAP simulations is required 
before one can attach cosmological significance to these 
findings. 

6. CONCLUSIONS 

We have generalized a previously described CMB 
Gibbs sampler to allow for exact Bayesian analysis of 
any anisotropic universe models that predicts a sparse 
signal harmonic space covariance matrix. This general- 
ization involved incorporation of a sparse matrix library 
into the existing Gibbs sampling code called "Comman- 
der" , and implementation of a new sampling algorithm 
for the anisotropy parameters given a sky map, P{6\s). 

We then considered a s pecial case of an i sotrop ic uni- 
verse models, namely the lAckerman et al.l (|2007t ) model 
which generalizes the primordial power spectrum P(fc) 
to include a dependence on direction, P(k). Explicit ex- 
pressions for the resulting covariance matrix is provided 
in their paper. 

We implemented support for this model in our codes, 
and demonstrated and validated the new tools with ap- 
propriate simulations. First, we compared the results 
from the Gibbs sampler with brute-force likelihood eval- 
uations, and then verified that the input parameters were 
faithfully reproduced in realistic WMAP simulations. 

Finally, we analyzed the five-year WMAP temperature 
sky maps, and presented for the first time the WMAP 
posteriors of the ACW model. The results from this anal- 
ysis are highly intriguing, but we emphasize that the ef- 
fect of instrumental systematics, particularly in the form 
of correlated noise, must be better understood before the 
findings can be given a cosmologically interpretation. 

Taken at face value, we find a preferred direction in 
the W-band WMAP temperature data pointing towards 
(/,6) = (110°, 10°) (Galactic longitude and latitude), 
with an anisotropy amplitude of 5* = 0.15 ± 0.039, for- 
mally corresponding to a 5.8a detection of 5, 7^ 0. Simi- 
lar results for 5, are found for the V-band data, although 
with a somewhat lower significance (g* = 0.10 ± 0.04; 
2.5tT). The preferred direction is very stable with re- 
spect to both data set and multipolc range. Figure [T2| 
illustrates the underlying anisotropic contribution for a 
simulation with parameters corresponding to the W band 
posterior. 

We have not been able to identify a plausible expla- 
nation for this effect in terms of known systematics. 
First, foregrounds do not appear to have much impact 
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Fig. 10. — Marginal ACW posteriors obtained from the non-template correeted V-band, P{h\d) (right) and P(gt\d) (left). Notice how 
P((/*|d) is shifted insignificantly with respect to the template-corrected V-band posterior. 
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Fig. 11. — Top: Posteriors for obtained from two 1-year WMAP W4 simulations with correlated noise. Bottom: The preferred 
direction posterior for one of the two above simulations. The peak positions of the second simulations are indicated by red dots, marked 
by (1) and (2). 



on the results, as consistent results are obtained both 
from foreground-corrected and raw maps. Second, al- 
though correlated noise does lead to a signature similar 
to the ACW model, its amplitude appears too low in 
the 5-year data. The least well constrained possibility 
is that of asymmetric beams, for which we lack proper 
simulations. 

While this particular signature certainly is highly in- 
triguing, we would like to point out that the main pur- 
pose of this paper is the demonstration of a general 
framework for analyzing anisotropic signal models. This 
is useful both for studying particular universe models 
(e.g., the ACW model), but also for understanding sys- 
tematic effects (e.g., correlated noise) in a given data set. 
We therefore believe that these methods may be useful 
in a wide range of applications, only some of which have 



been demonstrated in this paper. 
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red dots. Notice the rotational structure about the preferred direction. The amplitude of the anisotropic component is ~ ±15/iK, or ~ 3% 
of the isotropic component. 
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